function [pval,WaldBoot,Wald] = IPCA_gamma_delta_pvalue(X,W,Nt,Gamma,F,PSF,NBoot)
% Testing for significance of GammaDelta

% X is L*T
% W is L*L*T
% Gamma is L*(K+M)
% F is K*T
% PSF is M*T

df = 5;
L = size(X,1);
[K, T] = size(F);
GammaBeta = Gamma(:,1:K);
GammaDelta = Gamma(:,K+1:end);
M = size(PSF,1);

d = NaN(L,T);
FitDelta = NaN(L,T);
FitBeta = NaN(L,T);
for t = 1:T-1
    Wt = W(:,:,t);
    Ft1 = F(:,t+1);
    PSFt1 = PSF(:,t+1);
    Xt1 = X(:,t+1);
    FitDelta(:,t+1) = Wt*GammaDelta*PSFt1; 
    FitBeta(:,t+1) = Wt*GammaBeta*Ft1;
    d(:,t+1) = Xt1 - FitDelta(:,t+1) - FitBeta(:,t+1);
end

Gamma_Old = [GammaBeta zeros(L,M)]; F_Old = F;
WaldBoot = NaN(NBoot,1);

parfor b = 1:NBoot
    TimeIndex = randsample(2:T,T-1,true)';
    dBoot = [NaN(L,1) d(:,TimeIndex)];

    Student = sqrt((df-2)/df)*trnd(df,[1 T]); % unit variance and 5 degrees of freedom
    dBoot = bsxfun(@times, dBoot , Student);

    Xboot = FitBeta + dBoot;

    GammaBoot = IPCA_internal(Xboot,W,Nt,PSF,Gamma_Old,F_Old);
    GammaDeltaBoot = GammaBoot(:,K+1:end); GammaDeltaBoot = GammaDeltaBoot(:);
    WaldBoot(b) = GammaDeltaBoot'*GammaDeltaBoot;
end

GammaDelta = GammaDelta(:);
Wald = GammaDelta'*GammaDelta;
pval = sum(WaldBoot>Wald)/NBoot;
return

%--------------------------------------------------------------------------
function Gamma_New = IPCA_internal(Xboot,W,Nt,PSF,Gamma_Old,F_Old)
MaxIterations = 1000;
Tolerance     = 1e-6;
tol           = 1;
iter          = 0;
while iter<=MaxIterations && tol>Tolerance
    [Gamma_New, F_New] = IPCA_base(Gamma_Old,Xboot,W,Nt,PSF);
    tol       = max([ abs(Gamma_New(:)-Gamma_Old(:)) ; abs(F_New(:)-F_Old(:)) ]);
    F_Old     = F_New;
    Gamma_Old = Gamma_New;
    iter      = iter+1;
end
return